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1 INTRODUCTION 

• I— I 

Cold Dark Matter (CDM) is thought to dominate struc- 
5_j ture formation through gravitational clustering. As its name 
implies, this material starts out very cold, so that to an 
excellent approximation particles are confined to a three- 
dimensional surface [r, p(r)] in six-dimensional phase space. 
The distribution function / is infinite on this surface and 
zero elsewhere. The fine-grained distribution function is con- 
stant because it obeys the coUisionless Boltzmann equation. 
Hence, very high phase-space densities are in principle pos- 
sible right up to the present epoch. 

Studies of structure formation rely heavily on A'^-body 
simulations to follow the clustering of dark matter (DM). 
Galaxies are believed to form from baryons that have fallen 
in to the centres of DM haloes, so the implications of various 
theories of DM for the observable properties of galaxies de- 
pend sensitively on the small-scale structure of dark haloes. 
Hence, it is important to understand how accurately A'^-body 
simulations represent the inner parts of dark haloes, and this 
problem has attracted considerable attention (Jing et al. 
1995; Klypin et al. 2001; Taylor & Navarro 2001; Navarro 
2003; Power et al. 2003; Fukushige, Kawai & Makino 2003; 
Hayashi et al. 2003). 

On sufficiently small scales it is clear that the simula- 
tions will become untrustworthy because they employ only 
a finite number of particles (Kuhlman et al., 1996; Melott 
et al., 1997; Splinter et al., 1998; Binney & Knebe, 2002). 



ABSTRACT 

An estimate of the convergence radius of a simulated CDM halo is obtained under 
the assumption that the peak phase-space density in the system is set by discreteness 
effects that operate prior to relaxation. The predicted convergence radii are approxi- 
mately a factor 2 larger than those estimated for numerical convergence studies. 

A toy model is used to study the formation of sheets of the cosmic web, from 
which DM haloes form later. This model demonstrates the interplay between phase 
mixing and violent relaxation that must also be characteristic of spherical collapse. 
In the limit that sheets contain arbitrarily many particles, it seems that power-law 
profiles are established in both distance and energy. When only a finite number of 
particles is employed, relaxation is prematurely terminated and the power laws are 
broken. In a given simulation, the sheets with the highest peak phase-space densities 
are those that form from the longest waves. Hence simulations with little small-scale 
power are expected to form the cuspicst haloes. 
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Generally this scale is determined by resimulating a given 
halo with improved resolution, in terms of the number of 
particles employed, the number of waves used to define the 
initial conditions, the smallness of the gravitational soften- 
ing parameter e, and the smallness of the timesteps. Based 
on such simulation series. Power et al. (2003) define a ra- 
dius Tconv outside which a given simulated halo should be 
reliable, by calculating the smallest radius at which three 
criteria are all met: (i) the local orbital time is much longer 
than the numerical timestep; (ii) the centripetal acceleration 
of a circular orbit is smaller than O.Sl^oo/^i where V200 is 
the halo's characteristic velocity; and (iii) the local two-body 
relaxation time is longer than the Hubble time. 

While it is clear that satisfaction of all three criteria is a 
necessary conditions for a N-body simulation to be credible, 
it is not clear that it is a sufficient condition: errors in the 
integration of the equations of motion at any cosmic epoch 
can invalidate the final halo. We argue here that a particu- 
larly important epoch to get right is the one that ends when 
overdensities first collapse to form sheets that are virialized 
in only one direction. We present evidence that discrete- 
ness effects at this primary stage, rather than, for example, 
two-body relaxation in virialized haloes, are what limit the 
central densities of simulated haloes. 

We use a toy model and analytic calculations to clarify 
what happens when virialized CDM structures form, and to 
understand the extent to which N-body simulations of struc- 
ture formation are compromised by artifacts due to discrete- 
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ness. The model explains why cuspy profiles form even when 
power at high spatial frequencies is removed from the initial 
fluctuation spectrum. It also suggests that in the absence 
of discreteness effects, virialized haloes would have power- 
law cusps. Finally analytic arguments are used to estimate 
the radius in a simulated halo within which the central pro- 
file is flattened by discreteness effects that operate prior to 
virialization, and we show that these estimates are in good 
agreement with the latest numerical convergence studies. 



2 A TOY MODEL 

In this section we use a one-dimensional toy system to study 
the interplay between phase mixing and violent relaxation, 
and to show how discreteness cuts short this interplay. The 
simulations of Melott (1983) are similar, except that they 
are for hot dark matter and they employ a conventional 
particle-mesh N-body technique to advance the particles. 

Consider the dynamics of an odd number 2n -|- 1 of flat 
mass sheets that are all perpendicular to the x axis and move 
parallel to this axis. To avoid confusion later, we henceforth 
call these sheets 'particles'. The particles, which each have 
surface density E/n, interact only through gravity. There 
is no net force on the middle particle. Let k be the rank 
of a particle with respect to this middle particle, with k < 
for particles on its left. Then the particle with rank k 
experiences a force Fk — ~4nG'E{k/n) that only changes 
when the particle passes through another particle. Hence 
the particle's location, Xk{t), is quadratic in time between 
such passages. So we can analytically solve for the instant at 
which the next passage occurs, and update the phase-space 
coordinates of the particles to this time to machine precision. 
The ranks of the passing particles are then exchanged and we 
solve for the instant of the next passage. Hence this model 
provides a convenient laboratory for a high-precision study 
of virialization. 

We adopt initial conditions appropriate to the moment 
of turnaround of an overdense region of the Universe. The 
ranking of the particles will not have changed up to this 
time, so the requirement that at some fixed time r in the 
past, all particles were located at the origin becomes 



a;fc - UfcT - 27rGE(fc/n)r^ = (fc = 



(1) 



This set of equations has a solution that corresponds to a 
homogeneous distribution at turnaround 

a;fc = 27rGSr^(fc/n), Ufc = 0. (2) 

If we define = Xk ~ x^-i — /n and Vk = Vk — 

Vk-i, then equations (0 yield a relation between Xt and 
Vk, namely 

Vk^^. 



(3) 



FiglTlshows from top left to bottom right six stages in the re- 
laxation of 301 particles from the initial conditions obtained 
by choosing 



Xk = -0.6— cos(0.37rfc/n). 



(4) 



The time, in units of r, since these initial conditions were 
imposed is given in the top left corner of each panel. The 
negative slope in the top left panel implies that at the start of 
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Figure 2. The number of particles with \x\ less than a given value 
at three times. The bold lino shows the slope for N oc 



the simulation, the system is already collapsing. In the next 
panel {t = 0.4) the vertical orientation of the distribution 
near the origin indicates that the centre has at that time 
finished collapsing. By the time, t — 1.64, of the next panel, 
the centre has expanded and collapsed once more. The edge, 
by contrast, is still in full collapse. The following panels show 
that the phase lag ip between the oscillations of the centre 
and the outside grows rapidly. By the time, t = 9.7, of the 
last panel, is many tt, and near the centre it has become 
hard to follow the spiral of phase points. 

Fig. 121 explores the evolution of the density profile of 
this system by plotting the cumulative number of particles 
with \x\ less than any given number at three times. In this 
logarithmic plot, the curve for t — has unit slope to within 
the errors because the density of particles p ~ constant ini- 
tially. The slopes of the other two curves are near i, which 
implies p ~ x~^^^. 

By Jeans' theorem, the distribution function of the 
equilibriated system must be a function f{E) of energy 
E = ^v^ + $. From the fact that the force law is 



F{x) = -47rG'E(x/a 



\l/2 



(5) 



where a;niax is of order the amplitude of oscillation of the 
least bound particles, it is easy to show that the power-law 
nature of the density profile implies that 



f 



where 



/o = 2 



-3/2 



\ 3r2 



1/3 



dy 



(y2 + 1)5/6 



(6) 



(7) 



Two processes are manifest in Fig. (i) The winding 
up of the line of phase points into an ever tighter spiral is 
phase mixing, (ii) Violent relaxation in the sense of Lynden- 
Bell (1967) causes the edge of the occupied part of the phase 
plane to move outwards, and the points near the centre to 
move towards the origin: the system's time-varying gravi- 
tational field is transferring energy from the central to the 
peripheral particles. The process by which this energy trans- 
fer occurs is this. Between the top left and top right panels, 
the particles that are at \px \ ^0.2 in the second figure have 
fallen together under their mutual self gravity. Particles fur- 
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Figure 1. From top left to bottom right, six stages in the violent relaxation of 301 particles. The number in the top left corner of each 
panel is the time since the initial conditions of equations were imposed. The units of time are the turnaround time t of the underlying 
homogeneous distribution of particles. 



ther out have not been affected by this infall. However, when 
the inner particles expand, the outer particles are falling in 
past them, and the inner particles have to climb out of a 
deeper well than they fell into. Conversely, the outer parti- 
cles fall into a well that has significantly weakened by the 
time they rise up the other side. 

This transfer of energy from the inner to the outer par- 
ticles increases the density contrast between the centre and 
the outside. Since the frequency of a particle's oscillations 
through the centre scales as the square root of the mean den- 
sity p interior to it, the energy transfer enhances the rate at 
which the function ^1>{E) is steepening. The energy transfer 
works most effectively between groups of particles that differ 
by ~ 7r/4 in tp. So, as the first five panels of Fig. clearly 



show, the characteristic distance scale of the transfer rapidly 
decreases. In a useful, if oversimplified^ picture, the first col- 
lapse transfers energy between the inner and the outer half 
of the particles. The second collapse at the centre transfers 
energy from the first quartile to the second quartile, and a 
little later energy is transferred from the second to the third 
quartile, and later still on out to the fourth quartile. By this 
time the additional central collapse pictured in the fourth 
panel of Fig.Qis transferring energy from the innermost | of 
the particles, and so on. At each stage the transfers become 
smaller because they involve a smaller and smaller fraction 

^ The actual numbers of particles (out of 301) losing energy at 
successive central collapses are ~ 100, 44, 28 and 19. 
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Figure 4. As Fig.^but zoomed in on the centre of phase space and at comparatively late times. 



of the mass, but Fig. |3 shows that they remain significant 
to remarkably late times. 

Fig. 0] shows zoomed-in views of the centre of phase 
space at six fairly late times. In the first three panels the 
spiral of phase points is unmistakable, but in the last three 
panels it has become hard to trace because it bends through 
a large angle between points. In these circumstances violent 
relaxation can no longer transfer energy between particles, 
and the increase in the central mass density ceases. Fig. |3] 
confirms this conclusion by showing that the energy of the 
innermost particles stops decreasing at about the time, t = 
4.04, of the bottom left panel of Fig.El 

The last panel of Fig. |1] gives the impression that the 
system has a finite central phase-space density. If we re-ran 
the simulation with more particles, we would be able to trace 



the spiral even at this time, and the impression of a finite 
central phase-space density would be dispelled. Violent re- 
laxation would continue to later times and generate a higher 
central mass density. Fig.j^shows that once the phase-space 
spiral has ceased to be resolved, the energy of the inner- 
most particles increases, presumably as a manifestation of 
two-body relaxation (i.e., stochastic heating). 

In is interesting to quantify the peak phase space den- 
sity. This will be of order F"'^, where F is the Poincare in- 
variant of the orbit of particles with rank |fe| ~ 2. From the 
force law Q the energy equation is 

ii^ = 4^GE— ^ (X^/^ - x'/^) , (8) 
where X is the amplitude of a particle's oscillations. Equat- 
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Figure 3. Energy of the innermost 8 particles (full curve) and 
outermost 8 particles (dotted curve) as a function of time. For 
clarity the energy of the inner particles has been multiplied by 
2000. 



ing r to the product of X and the corresponding peak speed 
for orbits with Ifcl = 2, we find 



/max p — {3 T^GYiXy^. 



\-i/2 (n\ 

^) [2) 



7/2 



(9) 



This equation can be exploited in two ways. First we can 
consider the effect of increasing the numerical resolution at 
which a given structure is simulated. Then n increases while 
the other factors are constant, so the peak value of the DF 
rises steeply with n. 

Alternatively, we can apply equation to a series of 
structures within a given simulation. Then E and a;max are 
both proportional to n, so overall /max ~ n'^^^. Hence the 
highest phase-space densities are to be found in the largest 
structures, because it is in these structures that violent re- 
laxation works most efficiently. 



2.1 Case of finite initial phase-space density 

Consider the case in which the initially populated region of 
phase space has a finite width - as in the case of warm dark 
matter. The winding up of the spiral and violent relaxation 
will proceed as before until the spiral's radius of curvature 
becomes comparable with the width of the populated sheet. 
At this point violent relaxation ceases to be possible. Hence, 
the peak phase-space density of the relaxed system will be 
comparable to the initial phase-space density of the sheet. 
Melott (1983) reached a similar conclusion on the basis of 
his simulations of pancake formation in a hot-dark-matter 
cosmology. 



2.2 Summary of experimental results 

From our one-dimensional example we conclude the follow- 
ing 

(i) Material at each distance from the centre executes os- 
cillations. The phase of these oscillations becomes a progres- 
sively steeper function ip{E) of distance from the centre of 
phase space. 



(ii) Gravity transfers energy with great efficiency between 
particles whose oscillations are out of phase by ~ 7r/4. Even 
a small initial enhancement of density towards the centre 
causes this transfer to be outward. 

(iii) The outward transfer of energy enhances the density 
contrast, which steepens the function ipiE), and reduces the 
length scale over which individual energy transfers occur. 

(iv) Outward energy transfer continues until the underly- 
ing spiral in phase space ceases to be resolved in the sense 
that the spiral turns through a large angle between phase 
points. 

(v) In the limit of an arbitrarily large number of particles, 
the spiral remains resolved through arbitrarily many wind- 
ings, and the density profile of the final relaxed structure is 
strictly singular. 

(vi) The large dynamic range between the smallest 
lengthscale of energy transfer and the overall system size, 
combined with the hierarchical nature of the energy trans- 
fer process serve to establish a power-law density profile p ~ 
x"^^^ and a coarse-grained phase-space density / ~ E~"'^^ . 

(vii) Once the phase-space spiral has ceased to be re- 
solved, or its width is no longer negligible, a core is estab- 
lished and two-body relaxation causes the halo to heat the 



3 RELEVANCE FOR COSMOLOGY 

Has the toy model any relevance for realistic three- 
dimensional systems? The interplay between phase mixing 
and violent relaxation that it so clearly illustrates is cer- 
tainly generic and as such relevant to nearly spherical sys- 
tems. The period at which a star oscillates radially in a 
spherical potential is strictly a function of both energy and 
angular momentum. However, in a typical galactic potential 
the radial period of a star is largely determined by the star's 
energy.^ Consequently, motion in the {r,Pr) phase plane is 
closely analogous to the motion in the {x,px) plane stud- 
ied in the last section. The shells or ripples that are seen 
in unsharp-masked images of a significant fraction of ellip- 
tical galaxies (Malin & Carter 1983) are manifestations of 
the winding up of an initially narrow distribution of stars in 
the {r,pr) phase plane (Hernquist & Quinn 1988) very much 
like what we see in Fig. Q 

As Zel'dovich (1970) pointed out, the non-linear stage 
in cosmic structure formation is initially a one-dimensional 
process. It gives rise to the network of very overdense sheets 
that are now familiar as the 'filaments' of the 'cosmic web' 
that are seen in slices cut through computer simulations of 
gravitational clustering (Shandarin et al., 1995). The forma- 
tion of such a virgin sheet will closely parallel the dynamics 
of the last section. Differences of detail will arise early on 
through the continued transverse expansion of the matter, 
which makes the force acting between two wafer-thin slices 
of matter an explicit function of time, but these differences 
are likely to be unimportant, given the finite duration of the 
one-dimensional relaxation process, and the tendency of the 
motion in the perpendicular directions to be faltering and 



^ The dependence on L vanishes entirely in the case of the 
isochrone potential. 
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turning to contraction soon after a singularity is reached in 
the first direction. 

Thus virgin sheets of the cosmic web will have singular 
central densities p ~ a;~^/^. The characteristic transverse 
length scale of these first sheets depends upon the power 
spectrum of the initial fluctuations. Unless the power spec- 
trum has a sharp cutoff, many different length scales will 
characterize the transverse structure of the virgin sheets 
within a sufficiently large volume of space because many 
different linear modes combine to determine how and when 
a given part of the Universe will go non-linear and collapse. 

Eventually the sheets collapse transversely, either upon 
themselves or through falling onto the massive objects that 
tend to form at the nodes of the original cosmic web. During 
this second stage of structure formation, violent relaxation 
can occur again. It may be possible to idealize the dynamics 
as occurring in the phase plane defined by the direction of 
fastest infall, y, and its conjugate momentum. In this plane 
an initially narrow distribution of phase points will be wound 
up into a spiral, and the final density profile is likely to 
be a cusped function py ~ with /3 > 0. The strong 
initial inhomogeneity of the system will result in the spiral 
winding up much more rapidly than it does when a virgin 
sheet forms, so we expect /3 < | . Finally the over-dense tube 
resulting from collapse in the x and y directions will collapse 
along its length, z. Again a cusped final density profile in z 
is to be expected, pz ~ z but violent relaxation will not 
be efficient, so 7 < i. The final density profile p(r) can be 
thought of as the spherically averaged product of the profiles 
established in each collapse direction. It is easy to show that 
p j--(°'+f>+t) ^ where a = i. The exponent of r must thus 
be greater than ^ and is expected to be significantly less 
that |. Simulations of cosmic structure formation favour 
exponents that lie towards the upper half of this expected 
range (Hayashi et al. 2003). 

3.1 Discreteness in cosmic N-body simulations 

We have argued that when a sheet of zero thickness in phase 
space collapses, the final virialized structure should have a 
singular central density. We have also shown that there is a 
maximum phase-space density that can be achieved when a 
finite number of particles is used to sample the phase-space 
sheet [eq. ||5J]. Here we ask what the maximum phase-space 
density is that can be achieved in the best current N-body 
simulations of cosmic structure formation, and inside what 
radius within a DM halo this restriction becomes important. 

We have seen that one of two effects can set an upper 
bound on the phase-space density that occurs in a simula- 
tion. One is due to the termination of violent relaxation by 
discreteness. For the one-dimensional case it is quantified 
by eq. Q . The other is set by any upper limit on the initial 
phase-space density of the initial conditions, since Liouville's 
theorem forbids any increase in the phase-space density. 

It is not evident how to extend eq. @ to the three- 
dimensional case. This is unfortunate, for its implication 
that the highest phase-space densities are associated with 
the most massive objects, is unexpected. It also needs clari- 
fication: it does not apply to objects that form through the 
merger of virialized pieces. Hence, when there is considerable 
small-scale power, few if any objects to which it applies will 
have survived to the present epoch. However, the result does 



explain why Knebe et al. (2002) found that suppression of 
small-scale power in the input spectrum did not reduce the 
cuspiness of the final haloes. Also, it suggests that the mate- 
rial that now forms the cores of simulated haloes virialized 
in narrow ranges of length-scale and redshift: this material 
was in the largest structures to collapse before substructure 
could form within them. 

In view of the difficulty of generalizing equation @ to 
three dimension, we concentrate on ways in which discrete- 
ness effects limit the phase-space density achievable prior 
to virialization. Below we show that the upper limit on / 
that we derive from these effects predicts convergence radii 
for DM haloes in simulations that are in reasonable agree- 
ment with those determined empirically. This finding sug- 
gests that restrictions on / prior to virialization are domi- 
nant in practice. 

Prior to virialization discreteness limits / in two ways: 

(i) With a small number of sampling points, it is not pos- 
sible accurately to trace a curved phase-space sheet. 

(ii) In simulations the gravitational forces are estimated 
by concentrating mass into lumps. This concentration adds 
high-frequency noise to the force field. Neighbouring parti- 
cles sample this noise in largely uncorrelated ways, and thus 
acquire uncorrelated components of momentum that are dis- 
creteness artifacts. These momentum components broaden 
an initially perfectly thin sheet of phase points. 

We estimate these effects in Appendix A. The second 
effect is estimated in two complementary ways. In every case 
we estimate the Poincare invariant F over which two or three 
particles are distributed in one pair of canonically conjugate 
coordinates. Since there are three such pairs, our estimate 
of the DF is 



where m is the mass of a simulation particle. 

Equations 1A5II . IIAf 4ll . and 1AI9II give three estimates 
of r. The expressions all contain the dimensional factor 
Ho{Aq)^, where Ho is the Hubble constant and Aq is the 
comoving separation of particles at early times. This factor 
is multiplied by a numerical factor that is comparable to 
unity, but probably a little smaller. Thus below we use the 
estimate 



' [//o(Ag)2]3' ^ ' 

where is a factor a little smaller than unity. Equation lA4t 
and equation ijAf 6|l independently imply that / ~ a-«/2 in 
the run up to virialization. 

When the system virializes, we have argued that the 
peak phase-space density will decrease only if the phase- 
space sheet is thinly sampled compared to its width. In Ap- 
pendix B we show that on virialization the peak value of the 
DF does not diminish much below that given by equation 
TTH . 



4 APPLICATION TO SIMULATIONS 

Hayashi et al. (2003) studied the convergence of the central 
structure of DM haloes in state-of-the art simulations of 
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structure formation. It is interesting to compare the conver- 
gence radius rconv that they derived empirically from these 
simulations with the one that follows our considerations. 

Let us assume that DM haloes have cores in which 
p{r) — po{r /ro)~^ and the velocity dispersion tensor is 
isotropic - simulations suggest that these assumptions are 
not seriously violated. Then, with the gravitational poten- 
tial taken to be zero at the centre, the distribution function 
of the system is easily shown to be 



f{E) = 



(12) 



(2S)5/2' 

We find the radius r-mi-n at which the DF from this formula 
at V = Q equals the maximum phase density of equation 
O- We find 

5/2 _ moi^iff 

327r5/2G3/2(por-o)i/277m' ' 

As one proceeds inwards from rmin, more and more of the 
density is contributed by parts of phase space at which the 
DF exceeds the value of eq. Ijllfl . In a simulation these por- 
tions of phase space are under-populated, so the density falls 
more and more below its true power-law value. Hence r-mi-n 
is our prediction of the convergence radius of a simulated 
halo. 

If A'^'^ particles are used to represent the matter in a box 
of side L, the mass of an individual particle is 

and Tniin satisfies 



(14) 



5/2 



47r3/2(Gporo)i''2Ar3^n^- (^^^ 

We apply this formula to the Gl sequence of models listed 
in Table 1 of Hayashi et al. For the relevant halo po'"o = 
6.34 X IQ'h-^ Mq(/i~^ kpc)"^ while L = 5/i"^ Mpc, so 

r^in = 2.14(ryf7,„/0.3)"^/'^(iV/256)"®''^ft"' kpc. (16) 

As one proceeds from A = 32 to A = 256 down the Gl 
sequence of models, we find that = 1.14, 

1.47, 1.53, and 1.53. So rmin is larger than rconv by a fac- 
tor ~ 1.5?7~^^^ that is only slightly greater than uinty. This 
agreement between our analytic estimates and the numeri- 
cal studies suggests that the peak densities within simulated 
haloes are indeed set by discreteness effects that operate be- 
fore virialization. 



5 CONCLUSIONS 

Semi-analytic work on the formation of cosmic struc- 
ture usually focuses on the collapse of spherical over- 
densities (Gunn & Gott 1972; FiUmore & Goldreich 1984; 
Bertschinger 1985; Subramanian, Cen & Ostriker, 2000). 
Actually one direction collapses before the other two, so 
the crucial first phase of the virialization process is one- 
dimensional (Zel'dovich, 1970; Shandarin et al., 1995). Fill- 
more & Goldreich (1984) obtained a self-similar solution 
of the collapse problem in this planar symmetry. Unfortu- 
nately, their results are not directly comparable to ours be- 
cause (i) they assumed that the transverse expansion con- 
tinued unperturbed, and (ii) their initial overdensity was 



constrained to diverge at the centre of the sheet. Conse- 
quently, for alH > their solution has a collapsed core, and 
the phase-space spiral has a time-invariant form: its scale 
changes with time, but it does not wind up in the manner 
of Fig. El 

Our toy model of the planar collapse of CDM reveals 
a strong interplay between phase mixing and violent relax- 
ation, in which an initially small phase lag of the oscillations 
between the centre and periphery of an overdensity causes 
energy to fiow outwards, and this outward energy flow has- 
tens the growth of the phase difference, which in turn causes 
the exchange of energy to repeat on ever smaller scales. The 
numerical evidence suggests that in the limit of arbitrar- 
ily many particles, power-law profiles in space and energy 
would be established as a natural consequence of the large 
dynamic range between the smallest scale of energy inter- 
change and the overall size of the system. Power laws are 
also consistent with the original infinite phase-space density 
persisting at the centre of the relaxed system, while admix- 
ture of 'air' ensures that the phase-space density is finite 
and outwards-decreasing elsewhere. 

When the system is represented by only a finite number 
of particles, random motions develop prior to virialization. 
Consequently, as collapse begins, the phase-space density is 
finite and the interplay between phase mixing and violent 
relaxation is prematurely terminated. The central phase- 
space density of the relaxed system is comparable to the 
phase-space density just before collapse, while further out it 
differs little from the value that would hold in the limit of an 
arbitrary number of particles. Consequently, the real-space 
density profile must turn over near the centre from the slope 
it would have in that limit. 

An interesting result from the one-dimensional model 
is that, in a given simulation, the peak phase-space density 
in simulated sheets of the cosmic web scales with particle 
number as / ~ n'^^^ . Hence, the sheets with the highest 
surface density achieve the largest phase-space density. The 
surface density of the sheets that form is increased if we sup- 
press small-scale power in the initial fluctuation spectrum, so 
there is a suggestion that simulations that have little small- 
scale power will produce the cuspiest haloes. This finding is 
consistent with the counter-intuitive results of Knebe et al. 
(2002). 

The extension to the formation of DM haloes of these 
results for the formation of the sheets of the cosmic web is 
not straightforward. Some insight can be gained by imag- 
ining that the original process repeats as each of the other 
directions experiences collapse. This idealization is problem- 
atic because the collapses of these directions are rarely well 
separated in time. However, it is probably legitimate to ar- 
gue that in the limit of infinite particle number, the final 
density profile should be a power law'' because the collaps- 
ing sheet brings no relevant scale into the relaxation process, 
and the latter has no intrinsic scale. Also it seems likely that 
violent relaxation will be less efficient in later collapses than 
it was in the first, because the phase spirals associated with 
these directions will wind up very quickly. This reasoning 



^ Figure 4 of Hayashi et al. (2003) argues against this conclusion, 
however. 
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allows us to set an upper limit | on the power-law index of 
the final density cusp, which is clearly greater than i. 

If we assume that in the limit of infinite particle number 
the slope would be a = 1, as in the NFW profile, then we 
can estimate the radius rmin at which the maximum phase- 
space density sampled equals the phase-space density that 
developed prior to virialization. This radius forms a reason- 
able estimate of the radius at which the simulated density 
profile turns away from a line of unit slope. We find that 
fmin is larger by a factor ~ 2 than the convergence radius 
Tconv of Hayashi et al. (2003) . This finding strongly suggests 
that processes that take place prior to virialization deter- 
mine the smallest radius at which a simulated halo can be 
trusted. 

The importance of knowing which structures in cosmo- 
logical simulations are reliable, is clearly very great. Brute- 
force convergence studies such as those of Power et al. (2003) 
do not constitute a straightforward tool for validating N- 
body simulations because the limit in which reality is ap- 
proached involves sending three quantities to zero: the mass 
of an individual particle, the softening length, and the wave- 
length of the highest-frequency mode included in the initial 
conditions. Before an entirely convincing convergence study 
can be made, one needs to know along which path in the 
three-dimensional space of these parameters one should ap- 
proach the origin (e.g.. Splinter et al. 1998). Hitherto there 
has been a tendendency for N-body simulators to focus on 
the impact that discreteness has on the dynamics of relaxed 
systems. This study suggests that dicreteness may have its 
dominant impact prior to virialization. 
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Figure Al. An underlying ZeI'dovich wave. 
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APPENDIX A: WIDTH OF PHASE SHEET 

The initial conditions for cosmological simulations are in- 
variably set up by superposing a number of plane-wave per- 
turbations to the density and velocity fields. Comoving co- 
ordinates are generally employed rather than the inertial 
coordinates employed in Section |21 In canonically conjugate 
comoving coordinates, Zel'dovich's (1970) analytic solution 
for the evolution of an isolated wave is 

x = q + F 

p = aHF, (Al) 

where H = a/a is the Hubble constant at scale factor a and 
we adopt 

F = ^sin(fcg). (A2) 

aok 

Here g is a Lagrangian coordinate and oq is the scale factor 
at which the wave breaks and a caustic forms. 

Al Effective width of a distorted sheet 

The phase-space density of a group of particles is calculated 
by evaluating the phase-space volume that they occupy. By 
definition, the latter is the product of three Poincare invari- 
ants Pi = J dxidpi, one for each index i, with pi the mo- 
mentum conjugate to Xi. Consider therefore Fig. lAll which 
shows the track in the {xi,pi) plane of a collapsing wave, 
and a number of particles that are intended to realize this 
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wave in a simulation. Any three consecutive particles define 
a triangle. We now evaluate the area of this triangle, which 
is the Poincare invariant that the particles define. 

Taking the modulus of the vector product {xi — X2,pi — 
P2) X (a;2— a;3,p2— Pa) of two differences between three vectors 
{x,p) of the form lAlll . we find the Poincare invariant 

r = ^[(gi - q2){F2 - F3) - (g2 - g3)(fi - F^)] (A3) 

If we now expand F\ and Fz as Taylor series around 52 , we 
find to leading order 

r= \Ha^F"{q2){q2 - qi){q-i - q2){qz - qi) 



2ao 



sin(fcq)(Ag)^ 



(A4) 



where in the last line we have assumed that the particles 
are uniformly spaced in q with interparticle separation Aq. 
In a flat universe, H = a~^''^Ho, where Ho is the Hubble 
constant at the current epoch, and we see from equation 
IIA4I I that r ~ a^/^. This estimate of F is credible up to 
the epoch of virialization, a — oq. Hence, at virialization the 
particles occupy a minimum Poincare invariant 



r = Ho{Aqf [i(feAg)ay''sin(fcg)]. 



(A5) 



This expression for the Poincare invariant is made up of a 
dimensional factor 



Fo = Ho{Aqf, 
and a dimensionless factor 
rj = i(feAg)ay^ sin(feg) 
c^^ikAqfa'J' 



(A6) 



(A7) 



where the second line follows because a wave breaks, and a 
high-density virialized structure forms, around g = 0. At the 
Nyquist frequency, fcAg = tt, so this factor will be of order, 
but less than unity is cases of interest. The other two factors 
are each ~ i. Hence, equation IIA6fl gives a useful order-of- 
magnitude estimate of the smallest achievable value of F. 



A2 Unequilibriated fluctuations 

The process of setting up a cosmological simulation can be 
conceptualized as a two-stage process. One first sets up a dis- 
crete representation of a perfectly homogeneous distribution 
of matter, and then one perturbs it by displacing particles 
in a systematic way that generally involves waves of some 
kind. Let us focus on potential fluctuations implicit in the 
unperturbed 'homogeneous' distribution of particles, since 
these fluctuations are clearly an artifact due to discreteness. 
The nature of these fluctuations is elucidated by the 'glass' 
approach to realizing a homogeneous matter distribution. In 
this technique (White, 1996) particles are first distributed 
randomly within the simulation volume and are then moved 
subject to the usual equations of motion with gravity made 
into a repulsive force. After a sufficient number of integra- 
tion steps, each particle comes to rest at the bottom of its 
nearest potential well. When gravity is restored to an attrac- 
tive force, the potential well becomes a hill with the particle 
at its summit. When particles are distributed on a lattice, 
they again arrive at the summits of hills, and the structure 
of these hills can be investigated by Ewald summation. 



When the 'homogeneous' configuration is perturbed, 
the hills are distorted, but only mildly if, as should be the 
case, the perturbation is imposed at sufficiently large red- 
shift for the imposed displacements to be small compared 
to the inter-particle separation. Since the structure of an in- 
dividual hill depends on the positions of all particles, each 
particle is displaced slightly from the summit of its own hill, 
and will roll down the hill as soon as the simulation is in- 
tegrated forward in time. Its rolling motion is nothing but 
a discreteness artifact and will cause the Poincare invari- 
ant associated with neighbouring particles to grow as PrAg, 
where pr is the momentum of a typical particle's roll. 

Let the parabolic form of the potential around the sum- 
mit be 

$0 2 



$(a;) = 



2{Aqy 



(A8) 



where <l>o is an appropriate constant. Then from the Hamil- 
tonian ^p^/a^ -|-<l>/o we easily find that x satisfies the equa- 
tion of motion 



x" + ^- 



$0 



(HoAqy 



0, 



(A9) 



where a prime denote differentiation with respect to a and 
an Einstein-de Sitter cosmology has again been assumed. 
This homogeneous equation has solutions a; cx a", where 



M 1+ / ^6^° 



(AlO) 



On dimensional grounds $o/(-ffoAg)^ ~ 1, so the value of n 
for the growing solution is 



1/2 



HoAq 



1 I HoAq 

4 'I" 1 /o "I" 



32$, 



1/2 



1/2 



HoAq 



(All) 



Inserting x — xo{a/ao)"'^ , where xo and ao are constants to 
be determined, into Pi = a^^^Hox' , we find that the Poincare 
invariant associated with neighbouring particles is 



F = xo^o [a/ao) ' . 



(A12) 



By Ewald summation we find that the potential l|A8^ 
governing the displacement of a single particle from its po- 
sition in the lattice has coefficient 



$0 = —{HoAqY 
2tt 



(A13) 



With this value of "I>o equations lAlll and IIA12II yield 



F = 



^ (a/ao)'''®'' (ffo Agxo) . 

ZTT 



(A14) 



If we set ao equal to the expansion factor at which virializa- 
tion occurs, then xo becomes the magnitude of the particle's 
displacement at this epoch, and will therefore be a signifi- 
cant fraction of Ag. Consequently, our new estimate of F is 
similar to that of equation 1A5II . 



A3 Graininess of the Potential 

Since the gravitational force is calculated from the positions 
of a finite number of particles rather than the density field 
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of the true continuum, it contains spurious high-frequency 
structure. A lower hmit on the magnitude of this structure 
is given by the force from the nearest neighbour - in a col- 
lisionless system, forces from neighbours should always be 
small compared to the overall force. The contribution that 
this makes to the spurious velocity of a particle is 



Sv 



dt 



GmAx 



a2[(Aa;)2 +e2]3/2 



(A15) 



where we have assumed that Plummer softening is being 
used and v is the true peculiar velocity. Substituting for Ax 
from equation HAIL and again assuming a flat cosmology, 
this becomes 
Gm 



Sv = 



Ho{AqY 



da [1 — (a/ao) cos(fcAg)] 



ai/2{[i _ (a/ao)cos(fcAg)]2 + (e/Ag)2}3/2 



(A16) 



The expression in braces on the bottom of the integrand 
starts near unity and tends to {e/ Aqf . Typically the soft- 
ening length e is a significant fraction of Ag, so the total 
variation of this term is not large. The square bracket on 
the top starts at unity and vanishes abruptly when a ~ ao. 
Hence, the dominant variation of the integrand comes from 
the factor on the bottom, and a reasonable approxima- 
tion is 



2Gma}/^ 

ov — 



(A17) 



Ho{Aqf ■ 

Multiplying by aoAx we obtain an estimate of the Poincare 
invariant generated by graininess of the potential 

2Gm 



r = {5v){aoAx) 



3/2r, 

-an 1 



(a/ao)cos(A;Ag)]. (A18) 



matter-dominated cosmology, Gm/Aq 



In a flat 

(//oAg)^/16, so this equation can be written 

3/2 

a/ao) cos(fcAg)], 



r = //o(Ag)^^[l 



(A19) 
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Figure Bl. Shaded region shows estimate of the Poincaree in- 
variant occupied by the central 5 points after virialization. 



makes the estimate inapplicable to the case of a densely 
sampled wave. We conclude that on virialization the peak 
phase-space density does not decline strongly. 



which is very again similar to equation lASt . 



APPENDIX B: EFFECT OF VIRIALIZATION 

Fig. IBll shows a phase-space distribution of points that is 
rather sparsely sampling a Zel'dovich wave as it collapses. 
From this point on the coarse-grained phase-space density 
will at most points decline rapidly as the sheet wraps round 
and round itself, trapping ever more air in the roll. So let 
us ask what the smallest Poincare invariant is that we could 
associate with the n particles nearest the centre of the roll. 
We take this to be the shaded rectangle. Its invariant is 

(jikAq)\ i/2sin(nfcAg) 



■Ho{Aq) 



k j " " k 

— ag' (kAq) 



(Bl) 



Again the dimensional factor Ho{Aq)^ emerges, but this 
time multiplied by a numerical factor that is larger by 
than the corresponding one in equation 1A5II . This factor is 
not large for a sparsely sampled wave, and violent relaxation 



